Spatio-temporal mapping of interictal epileptic discharges based on mutual information of concurrent EEG and fMRI

Introduction: Simultaneous acquisition of EEG and functional MRI (fMRI) data has proven to be a powerful technique to investigate the location and dynamics of epileptiform networks associated with interictal epileptic discharges (IED). However, in most analyses for IED, EEG dominates fMRI, i.e. the EEG determines the onsets of the IED events in the general linear model (GLM) fMRI analysis (1, 2). Furthermore, a particular shape of the haemodynamic response function (HRF) is usually assumed although the BOLD response to IED can significantly differ from the HRF in healthy controls (3, 4). This work demonstrates that an information theoretic analysis based on the mutual information (MI) between the EEG-score and the fMRI data provides further insight in the study of epilepsy with EEG-fMRI. The most prominent feature of MI is that it balances both imaging modalities. In addition, the approach does not require a-priori models for the HRF shape nor does it assume a linear relationship between the spiking activity detected on the scalp-EEG and BOLD responses (5). Our approach is further enhanced by nonparametric statistical assessment of the MI maps based on a 4D spatio-temporal wavelet packet resampling approach (6). The technique was evaluated in 5 epileptic patients, confirming neurological assessment in 4 cases and showing better or equivalent performance than GLM-based analysis. Methods: 5 epileptic patients were scanned with a Siemens 3T Trio MR scanner (approved by local research ethics comittee) with simultaneous EEG and resting state (eyes-closed) gradient echo EPI fMRI acquisition. MRI parameters (slices, TR, TE, flip angle, voxel size): BR with 3 runs (76, 31, 12 IED) with 400 images (20, 1.2s, 30ms, 90°, 3.75x3.75x6 mm3); LS with 1 run (462 IED) of 500 images, SB with 1 run (82 IED) of 1100 images, MW with 1 run of 1100 images (25, 1.5s,35ms, 85°, 3.75x3.75x5.5 mm3); CI with 1 run (280 IED) of 600 images (32 slices, 2s, 30ms, 90°, 3x3x3.75mm3); T1-w and/or T2-weighted images were also acquired for all subjects to facilitate anatomical localization. fMRI preprocessing: Datasets were corrected for motion and lowfrequency trends (high-pass filter with cut-off frequency of 0.008Hz) with AFNI, and voxel timeseries were z-normalized. EEG recordings-analysis: 32, 64 or 96 MR-compatible EEG cap (EasyCaps, FalkMinnow Services, Herrsching, Germany) was used according to the 10-20 system. Electrodes were equipped with an additional 5kΩ resistance and impedances were kept as low as possible. EEG was acquired at 5kHz using 1-3 BrainAmp MR compatible amplifiers (Brain Products, Munich, Germany) and was synchronized with the MR trigger. MR gradient and cardioballistic artefacts were removed from the EEG using Vision Analyzer (Brain Products, Munich, Germany) using average artifact subtraction methods (8). EEG data was subsequently downsampled to 250Hz, and IED were visually marked by an experienced electroencephalographer. MI analysis: For each run, we defined a timeseries (temporal resolution of 250Hz) with ON-OFF periods of epileptic activity using the onsets of the IED marked on the EEGscore and each IED was assumed of one TR duration. Each timeseries was then downsampled to the fMRI images (EEG-fMRI-score). The MI between the EEG-fMRI-score and the fMRI data was computed voxelwise based on the entropy of the fMRI timeseries itself and conditional on the EEG-score timecourse (5) at successive TR latencies (range between 0-15s) by shifting the EEG-fMRIscore timecourse. A kernel density estimator with uniform grid of 100 bins and bandwidth adapted to the number of images (7) was used to estimate the probability density functions. MI maps were thresholded using a nonparametric wavelet resampling approach. First, 19 surrogate datasets were created with identical autocorrelation structure to the original fMRI data (6) and analyzed using the same procedure. Second, the cumulative density function of the maximum statistic of the MI values under the null was estimated from the surrogate data and the MI threshold was computed (p-value < 0.0005) after multiple comparison correction across latencies (6). To condense the temporal information of the MI analysis, we created three maps with the maximum MI across latencies, a delay map with the latency corresponding to the maximum MI, and an HRF signal-change map with the amplitude of the HRF at the latency of maximum MI. As the computation of MI does not assume any particular shape for the HRF, this was voxelwise estimated via single-trial averaging of the BOLD responses to the IED. GLM analysis: For comparison, datasets were also analyzed with a GLM approach (3dREMLfit function in AFNI, three regressors created with the convolution of the IED events with the SPM canonical HRF and its temporal and dispersion derivatives, and the 6 translation and rotation realignment parameters). GLM maps were created displaying the t-value of the canonical HRF regressor for those voxels where the statistical significance of the full model (F-test) was equal to the one used to threshold the MI maps. Results and Discussion: Our results demonstrated for the first time that MI analysis considering the fMRI data and timing of IED can reveal focal BOLD changes concordant with electro-clinical localisation and, when available, with GLM-based findings, and allows exploration of the dynamics of haemodynamic changes at whole-brain scale. In this work, the localization of the IED focus mapped with MI analysis was concordant with electroclinical localization in 4 out of 5 patients (SB, LS, MW, CI). In one of these cases (patient SB) the MI analysis identified two clusters of activation in left middle frontal gyrus and superior frontal gyrus (Figure 1a,b), areas where the epileptogenic zone was presumed to be located. These clusters were not found with GLM analysis (Figure 1c) possibly due to deviations in the HRF model (Figure 1d), showing the great flexibility of this approach over standard GLM regarding the neurovascular coupling with no need for prior on the HRF. No conclusive results were obtained by either MI or GLM analysis for patient BR, possibly owing to a reduced number of IED. In datasets with more frequent IED (LS, MW and CI) we observed high concordance between both techniques and with electroclinical localization (Figure 2, case LS), supporting the use of a linear model for the BOLD response and the canonical HRF in these cases. In sum, MI analysis is a promising method for the study of epilesy with concurrent EEG-fMRI, especially in clinical cases where the HRF deviates strongly from its canonical form as in children or patients with brain lesions (9). References: (1) Vulliemoz et al., 2010. Epilepsia, 51(4):491-505; (2) Gotman, 2009. Epilepsia, 49(S3):42-51; (3) Bagshaw et al., 2005. Neuroimage, 24(4): 1099-112; (4) LeVan et al., 2010. Neuroimage, 50(1):15-26; (5) Fuhrmann-Alpert et al., 2007. Neuroimage, 34(4):1545-61; (6) Patel et al., 2006. Neuroimage, 31(3):1142-55; (7) Jones et al., 1996, J Amer Stat Ass, 91(433):401-07; (8) Allen et al., 2000. Neuroimage, 12(2):230-39; (9) Jacobs et al., 2008. Neuroimage, 40(2):601-14. Figure 1. Case SB: Maps of a) latency, b) HRF signal change for MI-analysis, and c) GLM maps, overlaid over the corresponding T2-w anatomical slice; d) HRF estimates in Anterior Cingulate Cortex (ACC), left Superior Frontal Gyrus (SFG) and left Middle Frontal Gyrus (MFG). MI-analysis reveals IED-related BOLD activations in SFS and MFG, not detected with GLM, concordant with the clinical workout.